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ABSTRACT 

A 2d unsteady bicharacteristic scheme with 
shock fitting is presented and its characteristic 
step, shock point step and boundary step are 
described. The bicharacteristic scheme is com- 
pared with an UNO scheme and the Moretti 
scheme. Its capabilities are illustrated by com- 
puting a converging, deformed shock wave. 

INTRODUCTION 

Commonly used difference schemes are capable 
of computing complex flow patterns and shock con- 
figurations, but they smear out shocks and other 
discontinuities. To avoid the effects of shock smearing, 
shock fitting has to be used. This is most naturally 
done by the method of bicharacteristics. 

Characteristic Normal Form 

Introducing the substantial derivate Dq = 
D t = dt + (wV), a directional derivative operator 
Db ■= (-BV) = d t + ((v — ag)V) along a bicharac- 
teristic, the abbreviations v g := vg, d g :— gX7 and 
the logarithmic pressure P and density g, the Euler 
equations can be written as 



-P>tP 

7 



D t P 
D t v = -^VP 



energy equation 
continuity equation 
momentum equation 



or, in characteristic normal form: 


characteristic differential equation 


direction 


D g =-D a P 

1 


(1) 


D B v g -^D B P =a(Vv-d g v g ) 


\ v-ag J 



where the characteristic differential equations (CDEs) 
are valid only on their specific characteristic surfaces, 
i.e. the particle path C — (1,^) and the Monge cone, 
generated by the bicharacteristics B — (l,v — ag). The 
unit vector g singles out a specific bicharacteristic on 
the Monge cone (fig. 1). 



Characteristic Step 

The bicharacteristic scheme was constructed on 
a Cartesian grid, in analogy to Hartree's method [1]. 
Using an educated guess of v and a in a grid point 
at the new time level the bicharacteristics are drawn 
backwards in the directions of the coordinate lines. At 
the intersection points of these bicharacteristics with 
the old time level, called 'footpoints', the quantities 
P, p, v are interpolated in second order from the grid 
points with a 2d nine point least-squares method (fig. 
1). Using the discretized CDEs, the quantities P, p, v 
then can be determined in first order at the new time 
level. A corrector step, which uses the first order values 
of v and a to reconstruct the bicharacteristics, yields the 
quantities P, p, v in second order, eliminating the lateral 
derivatives on the new time level according to Butler [2] . 




Figure 1: Characteristic Step 

Treatment Of Shocks 

To begin with, a shock point is moved using its 
previous velocity. At its new location, the (oblique) 
Rankinc-Hugoniot conditions are used to determine the 
downstream values. With these values, a bicharacteris- 
tic is drawn backwards antiparallel to the shock front 
normal (fig. 2). The correct position of the new shock 
point is found by iteration, such that both the Rankine- 
Hugoniot conditions and the differential equation along 
the bicharacteristic are fullfilled. 

Shock points are moved along their trajectories. 
The density of shock points on a shock front can be 
kept constant by creating or deleting shock points. Two 



its predecessor and successor. Though all shock points 
are stored in an arbitrary sequence in a single array, 
shock contours and their normal vectors can neverthe- 
less be reconstructed with the help of the pointers. New 
shock points are detected with an algorithm according 
to Moretti [3] at each new time level. 

Shock Point Step 



Shock Point 




Figure 2: Shock Point Step 



Boundary Step 

When a new time level is computed, it is first done 
regardless of all discontinuities. Then, the shock points 
are moved. Invalid grid points, which were computed 
by drawing bicharacteristics crossing a shock surface, 
are updated afterwards with the help of a boundary 
step: 

When a Mongc cone is intersected by a boundary, e.g. 
a shock wave, the bicharacteristics are no longer drawn 
backwards in coordinate directions but in directions 
normal and tangential to the boundary (fig. 3). This 
improves the accuracy of the step considerably. The 
intersection points of the bicharacteristics with the 
boundary are determined, the footpoint values are in- 
terpolated at the boundary surface, and the grid point 
is updated. The boundary step is only of first order 
because the lateral derivatives cannot be eliminated 
with the Butler procedure in case the bicharacteristics 
of a monge cone are cut off at different time levels. 




Figure 3: BoiinHa.rv Sten 



Interpolation 

In the vicinity of discontinuity surfaces the in- 
terpolation algorithm used for ordinary grid points 
becomes ill-conditioned, producing instabilities. The 
most preferable way of interpolation in this case proved 
to be a two-dimensional, second order least-squares 
pattern, sized 7x7 grid points, in which all upstream 
grid points are ignored, but in which all shock points 
of this area are taken into account. According to our 
experience an interpolation in smaller regions will 
render the scheme unstable. 

COMPARISON WITH OTHER NUMERICAL 
SCHEMES 

The bicharacteristic scheme was tested with an 
analytical solution, and computations of converging 
shock waves were compared with the corresponding 
results of a UNO scheme and the Moretti scheme (a 
A-scheme working with shock fitting). As expected, the 
bicharacteristic scheme turned out to be superior to the 
UNO scheme regarding the prediction of position and 
velocity of the shock front. The shock fitting algorithm 
provides exact information on position, direction of 
propagation, and local Mach number of the shock wave. 
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Figure 4: Comparision of convergence of the UNO-, A- 
and bicharacteristic scheme 

The Moretti scheme [4] proved to be comparable 
in accuracy to the bicharacteristic scheme (fig. 4). The 
aforementioned is faster (by a factor 1.5, assuming 
Courant Number 2.0). This is due to the time consuming 
interpolations in the bicharacteristic scheme. We found 
that the convergence of the shock front position towards 
a known solution by refinement of the grid is slightly 
faster in the bicharacteristic scheme, due to its bound- 
ary step. 

Recently, Nasuti and Onofri [5] extended Moretti's 
original shock fitting algorithm to handle triple points. 
We implemented their extensions in our version of the 



be published soon). We found, that the shock fitting 
algorithm proposed by these authors still has some 
disadvantages. For example: The fragmentation of the 
shock front, described in the next paragraph occurs 
too late and is not enough pronounced in case the 
shock contour has an unfavourable orientation to the 
computational grid; whereas our shock fitting algorithm 
is not influenced by the relative position of shock front 
and computational grid. 

COMPUTATION OF CONVERGING SHOCK 
WAVES 

The computation shown here (fig. 5 and 6) started 
with a slightly deformed shock wave of Mach No. 2.5, at 
radius 1.0. Isopycnics of a time step immediately before 
and some time after fragmentation occurs are shown in 
Figure 5 and 6, respectively. Presently, the computation 
of converging shock waves extends to the instant of 
fully developed fragmentation, just before reflection 
of the leading shock. The extension of the method of 
bicharacteristics to proceed beyond this point is under 
work and does not pose any fundamental problems. 

Calculating converging shock waves with a bichar- 
acteristic scheme, one has to make sure that the Mach 
stem will develop correctly. As mentioned above, infor- 
mation on the downstream values of the physical quanti- 
ties reaches a shock point by transportation along a ret- 
rograde bicharacteristic. The footpoint of this bicharac- 
teristic generally is located only fractions of a grid width 
apart from the shock front. However, the physical mech- 
anism which generates the Mach stem is a density hump 
which gradually steepens as the shock wave converges. 
Finally, this density hump takes the form of a bow shock 
[6, 7]. Emerging gradually from a compression wave, its 
shock profile is smeared out similar to shocks in common 
difference schemes over a number of grid points. For this 
reason the retrograde bicharacteristic sees only a slight 
increasing of e.g. density values where a marked den- 
sity hump should be located. Therefore, the local shock 
velocity could be calculated too small and, hence, the 
formation of the Mach stem could be delayed. 

To overcome this problem, the developing bow 
shock has to be detected with a pattern recognition 
algorithm, and consecutively be treated with the shock 
fitting algorithm. 

As mentioned above, the density hump gives rise 
to new developing shocks in the vicinity of the triple 
points (fig. 7). These lateral waves are not fitted yet 
because of the considerable effort, this might take. 
However, the computed solution is acceptable as long as 
the lateral waves remain sufficiently weak. Depending 
on our initialisation the Mach Number of the lateral 
waves was only slighly above 1.0 and therefore the 
results can be considered as valid. 
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Figure 5: Isopycnic at t=0.1750 (before fragmentation) 
The bold contour represents the leading shock wave. 
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Figure 6: Isopycnic at t=0.2622 (after fragmentation). 
The curvature of the Mach stem is below grid width 
and cannot be observed by shock capturing schemes in 
this resolution. Compare also fig. 7. 



OUTLOOK 

The following improvements seem feasible: 
-The treatment of intersecting shock point trajectories 
should be reconsidered. 

-Newly developing shock fronts and discontinuity lines 
should be detected and fitted as proposed by Moretti. 

The bicharacteristic scheme with shock fitting pre- 
sented here provides great accuracy and physical insight 
and allows a variety of applications: It could be used to 
compute, e.g., channel flows, flows around wings, MHD 
problems, star pulsation, and non-equilibrium flows. 
The scheme can also profitably be used as a standard 
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Figure 7: Density at t=0.2622 

Due to the shock fitting algorithm the leading shock 
wave, which is represented here by density values ex- 
actly fullfilling the Rankine-Hugoniot relations, is not 
smeared out. The eminent area between lateral shock 
wave and slip line consists of those fluid particles which 
passed both the leading shock wave and the lateral wave. 



NOMENCLATURE 



a 


sound speed 


B 


vector along a bicharacteristic 


Co 


vector along the particle path 


D t = d t + (vS7) 


substantial derivative 


D = d t + (vV) 


derivative along Co 


D B = (BV) 


derivative along B 


d g = (flV) 


derivative in direction of g 


d t = d/d t 


time derivative 


9 


unit vector (spatial) 


P = ln(p/p ) 


logarithmic pressure 


P,Po 


pressure, initial upstream press. 


V 


velocity field 


v g = vg 


projection of v on g 


7 


specific heat ratio 


Q = ln(p/p ) 


logarithmic density 


P, Po 


density, initial upstream density 
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